library(foreign)
library(readstata13)
library(stargazer)
library(ggplot2)

# Table A9
validation<- read.dta("Validation_crop.dta")
# State fixed effects
reg1 <- lm(lntotprod ~ monsoon + as.factor(stateno), data = validation)
summary(reg1)
reg2 <- lm(lntotprod ~ monsoon + monsoon2 + as.factor(stateno), data = validation)
summary(reg2)
# District fixed effects
reg3 <- lm(lntotprod ~ monsoon + as.factor(statedists), data = validation)
summary(reg3)
reg4 <- lm(lntotprod ~ monsoon + monsoon2 + as.factor(statedists), data = validation)
summary(reg4)

# Create table of results

omit.list <- c("stateno","statedists")
omit.labs <- c("State FE", "District FE")
omit.stats <- c("rsq","ser","f")
cov.labs <- c("Monsoon SPEI","(Monsoon SPEI)^2")

fit.list <- list(reg1, reg2, reg3, reg4)

sink("~/cropValidation.tex")
stargazer(title = "SPEI and Agricultural Production",
          header = F,
          covariate.labels = cov.labs,
          fit.list,
          omit.stat = omit.stats,
          omit = omit.list,
          omit.labels = omit.labs,
          no.space = T,
          suppress.errors = F,float = F,dep.var.labels = c("Total Agricultural Production"))
sink()

# Figure A4
# Generate a plot of ideal point for SPEI
p <- ggplot(data = data.frame(x = 0), mapping = aes(x = x))
fun.1 <- function(x) 10.63 - 0.48*x^2 + 0.67 * x
p <- p + stat_function(fun = fun.1) + xlim(-3,3)
p <- p + ylab("ln(Agricultural Production)") + xlab("Monsoon SPEI") + theme_bw()
p

